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A  FORTRAN  COMPUTER  PROGRAM  FOR  CALCULATING  THE  PROLATE 
SPHEROIDAL  ANGULAR  FUNCTIONS  OF  THE  FIRST  KIND 


INTRODUCTION 

The  Helmholtz  or  scalar  wave  equation,  for  steady  waves  (72  +  k2)>f>  *  0, 
where  k  »  2tt/X  with  X  equal  to  the  wavelength,  is  separable  in  prolate 
spheroidal  coordinates  (E,  n,  $).  The  factored  solution  is  written  as 
R(c,5)S(c,nH(<*>)  •  Here  c  *  ka/2,  where  a  is  the  interfocal  distance  of  the 
elliptic  cross  section  of  the  spheroid. 

The  angular  function  of  the  first  kind  S^^ (c,n)  is  one  of  two  indepen¬ 
dent  solutions  to  the  ordinary  differential  equation  in  the  angle  coordinate 
n  arising  from  the  separation  of  variables.  This  solution  is  characterized 
by  the  four  parameters  m,  £,  c,  n.  For  each  of  the  choices  of  m,  £,  c,  and  n 
there  exists  a  set  of  solutions  to  the  prolate  angular  equation,  each 
solution  characterized  by  a  separation  constant  or  eigenvalue  A^.  As  with 
the  corresponding  associated  Legendre  functions  of  spherical  geometry,  it  is 
often  convenient  to  specify  the  argument  n  in  terms  of  the  angle  0  *  cos  1 n ; 
i.e.,  S^^c.n)  =  S^£(c,cos0). 

The  computer  program  PANGFN  calculates  numerical  values  for  the  angular 
function  of  the  first  kind  S^p(c,n)  and  the  associated  eigenvalues  A_^  for 


m£ 


desired  values  of  m,  t,  c,  and  0.  PANGFN  is  intended  to  replace  the  prolate 
portion  of  the  FORTRAN  computer  program  ANGLFN  [1],  which  was  previously 
developed  at  the  Naval  Research  laboratory  (NRL)  to  evaluate  both  the  prolate 
and  the  oblate  angular  functions  of  the  first  kind.  Unfortunately,  ANGLFN 
and  two  companion  computer  programs  for  calculating  the  prolate  and  oblate 
radial  functions  of  both  kinds  and  tneir  first  derivatives  [2,3]  were 
developed  around  the  large  exponent  size  (±307)  of  the  CDC3800  computer  at 
NRL.  These  programs  are  not  easily  modified  to  run  on  computers  with  a 
significantly  smaller  exponent  range.  The  program  PANGFN,  however,  is 
designed  to  run  on  computers  with  any  exponent  range.  It  is  written  in 
universal  FORTRAN  and  should  run  on  any  computer  that  accepts  this  language. 
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Similar  universal  computer  programs  will  be  developed  in  the  future  for  the 
oblate  angular  functions,  prolate  radial  functions,  and  oblate  radial 
functions.  A  universal  program  called  LINPRO  [A]  already  exists  for 
calculating  the  linear  prolate  functions  and  eigenvalues.  These  functions, 
which  are  useful  in  the  representation  of  band-limited  and  time-limited 
physical  processes,  are  constructed  from  the  prolate  angular  functions  with 
m  set  equal  to  zero. 

ANALYSIS 

The  prolate  spheroidal  system  can  be  formed  by  rotating  the  two- 
dimensional  elliptic  coordinate  system,  consisting  of  ellipses  and  hyperbolas, 
about  the  major  axis  of  the  ellipse.  The  prolate  spheroidal  coordinates  (5, 
n,  with  1  i  5  -  ",  -l  in-  l,  and  0  1  $  S  2r  are  related  to  the  Cartesian 
coordinates  (x,y,z)  by  the  transformations: 


(a/2)[U2 

-  1)(1  -  n2)]j£cos$  , 

(1) 

(a/2)[U2 

-  1)  ( 1  -  n2)]^?^  , 

(2) 

and 


z  *  a£n/2  . 


(3) 


The  spheroidal  coordinates  5  *  const.,  n  *  const.,  and  4>  *  const,  define 
the  following  set  of  orthogonal  surfaces,  as  shown  in  Fig.  1: 


X^  4*  y^  2^ 

ellipsoid  of  revolution:  — 1 - 1 -  + - 


(a/2)2(C2-l)  (a/2)252 


(A) 


hyperboloid  of  two  sheets: 


x2  +  y2 


(a/2)2(l-n2)  (a/2)2n2 


-1  , 


(5) 


and 

half  plane  containing  the  z-axis:  d  *  tan  ^(y/x)  .  (6) 
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ris*  1  -  The  prolate  spheroidal  coordinate  system 


Separation  of  the  Helmholtz  equation 


(V2  +  k2)  '■}>  -  0 


in  prolate  spheroidal  coordinates  yields  the  following  solution: 

*u  ’  V(C|E)V(C',,)’.W  ■  (8> 

where  V<c,5)  is  the  radial  function,  S^Cc.n)  is  the  angular  function,  and 
$  ($)  is  the  azimuthal  function. 

m  n  f,.  f\  c  o'!  and  $  (6)  satisfy  the  following  ordinary 

The  functions  R^e.U,  s^c»n,»  aua  mw 

differential  equations: 

-  tv  - c^2  +  fitr'V " 0  •  <9> 


3 


(10) 


^ta-n2)-5^1  +  (Ami  -  cV  -  -  o  , 


—3  +  m2^,  „  o 

d$2  m  ’ 


(ID 


where  m  and  A  n  ara  the  two  separation  constants  or  eigenvalues  occuring  in 

the  separation  of  variables.  Throughout  the  rest  of  this  manuscript  the  term 

eigenvalue  is  understood  to  be  A^,  unless  otherwise  stated. 

In  physical  problems  in  which  the  field  has  to  be  periodic  and  unique 

over  the  range  of  the  azimuthal  coordinate  it  is  required  that  m  be  an 

integer.  For  fixed  m  and  c  +  0  the  numbers  A  ^  for  which  Eqs.  (9)  and  (10) 

have  non-trivial  convergent  solutions  are  ordered  numerically  in  an 

ascending  series  and  labeled  with  the  integers  £  *  m,  m+1,  m+2,....  When 

c — *K),  Eq.  (10)  reduces  to  the  ordinary  differential  equation  of  the 

associated  Legendre  functions  whose  eigenvalues  A^  are  equal  to  £(£+1). 

The  two  independent  solutions  of  Eq.  (9)  are  known  as  the  radial  function 

of  the  first  kind  R^P(c,0  and  the  radial  function  of  the  second  kind 

(2)  tiwl 

Rv  <>  (c,0 .  Similarly,  the  solutions  of  Eq.  (10)  are  known  as  the  angular 

m-t 

function  of  the  first  kind  S  p  (c,n)  and  the  angular  function  of  the  second 
(9)  m<- 

kind  S  o  (c,n).  Excellent  discussions  of  the  uses  and  properties  of  these 
me 

functions  are  found  in  the  monographs  by  Meixner  and  Schafke  [5]  and  Flammer 

[6]. 

Three  volumes  of  tables  of  numerical  values  for  the  prolate  radial 
functions  and  their  first  derivatives  with  resoect  to  5  were  published  in 
1970  [ 7 ] .  These  volumes  contain  entries  for  the  following  range  of  parameter: 
m  *  0  (Volume  1),  m  =  1  (Volume  2),  m  *  2  (Volume  3),  £  =*  m,  m+1,  ...,  m+49; 

5  *  1.00000001,  1.0000001,  ...,  1.01,  1.02  (0.02)  1.2*,  1.4  (0.2)  2.0,  4.0 
(2.0)  10.0;  c  *  0.1  (0.1)  1.0,  2.0  (1.0)  10.0,  12.0  (2.0)  30.0,  35.0,  40.0. 

A  single  volume  of  tables  of  numerical  values  for  the  prolate  angular 
functions  and  the  linear  prolate  eigenvalues  was  published  in  1975  (8).  The 
range  of  variables  covered  in  this  volume  is  m  =  0;  £  *  0  (1)  49;  9=0°  (1°) 
90°,  where  q  =  cosQ;  c  =  0.1  (0. 1)1.0,  2.0  (1.0)  10.0,  12.0  (2.C)  30.0,  35.0, 


The  notation  1.02  (0.02)  1.2  indicates  1.02,  1.04,  1.06, 
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40.0.  Three  volumes  of  oblate  radial  functions  [9]  and  one  volume  of  oblate 
angular  functions  [10]  were  also  published. 


PROLATE  ANGULAR  FUNCTION  OF  THE  FIRST  KIND 

Since  the  angular  function  of  the  first  kind  S^^CjH)  reduces  to  the 
associated  Legendre  function  of  the  first  kind  P^Cn)  in  the  limit  c—^O,  it 
is  convenient  to  expand  the  angular  function  in  the  following  series: 


The  expansion  coefficients  dn(c|m£)  are  sometimes  called  the  D  constants. 
The  prime  sign  on  the  sum  indicates  that  n  *  0,  2,  4,  ...,  if  £-m  is  even, 
or  n  *  1,  3,  5,  ...,  if  £-m  is  odd.  This  ensures  that  the  angular  function 
S^(c,n)  has  the  same  evenness  or  oddness  with  respect  to  n  as  the  correspond¬ 
ing  associated  Legendre  function  p£(n),  i.e.,  S^^(c,n)  an  even  function  of 
n  if  £-m  is  even  and  is  an  odd  function  if  £-m  is  odd.  Ferrers'  [11] 
definition  of  the  associated  Legendre  functions  has  been  adopted.  Thus, 


P®(n)  -  (1-n2) 


m/2 


This  leads  to  the  following  recursion  relation: 

(n-m+i)p“  , (n)  -  (2n+l)nP^(n)  -  (n+m)P^  , (n)  ,  (14) 

n+i  n  n-i 

where 

Pm  ,(n)  -  0  ,  (15) 

m-i 

Pm(n)  -  (2m-i) ! !  (l-n2)m^2  ,  (16) 

m  7 

with  (2m— 1) 1 !  a  (2m-l) (2m-3) . . . (3) (1) . 
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Substitution  of  Eq.  (12)  into  Eq.  (10)  and  use  of  the  above  recursion 
relation  for  the  Legendre  functions  leads  to  the  following  three-term 
recursion  relation  in  terms  of  the  D  constants. 


(2m+n+2)  (2m+n+l)  2  j 

(2m+2n+3)  (2ra+2n+5)  n+2 


+ 


|^(m+n)  (m+n+l)  - 


.  2(m+n)  (m+n+1)  -2nr-l  2I 

m£  (2m+2n+3) (2ra+2n-l)  C  J 


d 

n 


+  - B-fo-P, _  c2d 

(2m+2n-3) (2m+2n-l)  1 


0  . 


(17) 


Equation  (10)  has  two  regular  singularities  at  n  ■  ±1.  If  the  choice  of 
*4^,  is  arbitrary,  solutions  given  by  Eq.  (12)  may  be  divergent  at  either 
n  *  1  or  n  *  -1.  It  is  necessary  in  physical  applications,  however,  that 
S^^(c,n)  be  finite  at  both  n  *  l  and  n  *  -1.  It  is  also  required  that 
Eq.  (12)  converge.  This  requires  that  d  ,0/d  —►0  as  n — Use  of  eigen- 

riT\i  n 

values  that  satisfy  both  this  condition  and  Eq.  (17)  above  will  result  in 
successful  expansions  of  S^(c,n). 

The  desired  eigenvalues  are  obtained  by  a  two-step  process.  First  an 
approximation  to  the  eigenvalue  is  calculated  by  use  of  formulas  such  as  a 
polynomial  expansion  in  c.  With  this  approximation  as  a  starting  value, 
the  eigenvalue  is  then  calculated  using  a  variational  procedure  developed  by 
Bouwkamp  [12].  In  the  Bouwkamp  procedure  the  three-term  recursion  relation 
found  in  Eq.  (17)  is  rewritten  in  the  following  two  forms: 


N* 

n 


-e 


m 

n 


vm-A  n+N111 
‘  n  ra£  n+2 


,  n  >  2  , 


(18) 


and 


+  A  p  -  — 
me  „m 
N 

n 


,  n  >  2  , 


(19) 
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where 


m 


(m+n)  (m-hvfl)  +  (c2/2) 


2)  1- 


4m2 -1 


(2m+2n-l) (2m+2n+3) 


.  n  >  0  ,  (20) 


„m  n(n-l)  (2nrfn)  (2nvfn-l)c 
n  *  “(2m+2n-l)2  (2nrf-2n-3)  (2mf2n+l) 


♦  n  >  2  , 


(21) 


„m  (2m4n)  (2m+n-l) _  2  ^n 

n  *  (2m+2n-l) (2mr2n+l)  C  d  0 

n-4 


(22) 


Requiring  that  d^+7/dn — *"0  as  n — <-  *  is  equivalent  to  requiring  that  M01 — •*-0 

as  n— ►  «.  Equation  (18)* can  be  rewritten  as  a  continued  fraction  in  terms 

of  ym  »  •••»  Bn,  8m,7,  and  A  n.  Equation  (19)  can  also  be 

n  iiTi  •  n  n.'T'Z  tut- 

rewritten  as  a  continued  fraction  in  terms  of  y11  ,,  ym  ,,  S'11  Bm  , 

n-2’  n-4’  *  n-2’  n-4, 

. . . ,  and  using  the  fact  that  *  -y^  -5*  A^  and  »  -y™  +  A^  to 

terminate  the  sequence. 

The  starting  value  for  the  eigenvalue  A^  is  inserted  into  the  two 
continued  fractions,  one  with  diminishing  subscripts  and  one  with  increasing 
subscripts.  m+9  is  calculated  using  both  expressions,  and  the  difference 
in  the  two  values  is  used  to  determine  a  correction  to  The  process  is 

repeated  until  A^  is  obtained  to  the  desired  degree  of  accuracy. 

When  the  value  of  A^  has  been  accurately  determined,  the  D  constants 
can  be  calculated  by  successive  application  of  Eq.  (17).  The  D  constants 
obtained  from  Eq.  (17)  are  un-normalized  since  this  equation  is  homogeneous. 
A  suitable  normalisation  is  established  by  requiring  that  the  angular 
functions  have  the  same  normalization  as  the  associated  Legendre  functions. 
Thus 


"*  9  t"  2 

^ [sie)(c’n>]  dn  -  j[vt (n)]  dri  =  drlryli-'^!'  *  (23) 

Substituting  the  expansion  for  (c,p)  into  Eq.  (23)  and  using  the  known 
orthogonal  properties  of  the  Legendre  functions  provides  the  following 
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normalizing  relation  for  the  D  constants: 


00 


n-0,1 


(n+2m) !  d2 


n 


(2m+2n+l)  n! 


2  (£+m)  1 
(2-d+l)  (£-m)  l 


(24) 


This  normalization  scheme  was  first  used  by  Meixner  and  Schafke  [5].  It  has 
the  practical  advantage  of  eliminating  the  need  to  numerically  evaluate  the 
normalization  factor  J"  (c,n) ]2dn  which  is  often  encountered  in  problems 

involving  expansions  in  angular  functions.  Once  the  D  constants  have  been 
normalized,  the  angular  function  can  be  evaluated  using  Eq.  (12).  Alterna¬ 
tively,  the  angular  function  can  be  evaluated  using  un-normalized  D  constants 
and  then  corrected  by  multiplying  by  the  ratio  of  the  normalized  to 
un-normalized  value  of  any  one  of  the  D  constants,  say  d» 

■(•“HI 

BRIEF  DESCRIPTION  OF  PANGFN 

The  computer  program  PANGFN  calculates,  in  double  precision  arithmetic, 
numerical  values  for  A^  and  S^^(c,n)  for  desired  values  of  m,  &,  c,  and  6, 
where  6  *  cos-1(n).  The  program  is  written  in  universal  FORTRAN  and  should 
run  on  any  machine  that  accepts  the  language.  The  main  program  calls  four 
subroutines:  "PLEG",  "GETEIG",  "CONVER",  and  "OUTPUT”.  "PLEG"  generates  the 
associated  Legendre  functions;  "GETEIG"  produces  an  approximation  to  the 
eigenvalue  for  starting  the  Bouwkamp  procedure;  "CONVER"  uses  the  Bouwkamp 
procedure  to  determine  the  eigenvalue;  and  "OUTPUT"  prints  the  final 
results.  The  processes  of  these  routines  are  described  in  more  detail  in  the 
section  entitled  COMPUTATIONAL  PROCEDURES. 

PANGFN  is  designed  to  accommodate  the  exponent  range  and  word  length  of  the 
user's  computer.  It  is  necessary  to  change  the  first  two  executable 
statements  in  the  program  to  correspond  to  the  user's  computer.  These 
statements  specify  NDEC,  the  number  of  decimal  digits  available  for  double 
precision  numbers,  and  NEX,  the  maximum  possible  exponent.  In  the  unlikely 
event  that  NDEC  exceeds  36,  the  value  for  t  given  by  PI  in  the  third 
executable  statement  must  be  extended  to  NDEC  digits.  The  user  should  also 
set  the  array  dimensions  large  enough  to  accommodate  the  desired  parameter 
ranges,  as  described  in  the  section  entitled  DIMENSIONS  AND  STORAGE. 


8 


The  remainder  of  this  report  describes  the  program  PANGFN.  Included  are 
a  listing  of  the  significant  FORTRAN  variables  and  a  description  of  the  major 
computational  blocks.  A  discussion  of  parameter  input  and  resulting  output 
follows.  A  listing  of  the  program  and  a  sample  output  are  also  given. 


SIGNIFICANT  FORTRAN  VARIABLE  NAMES 

ARG:  Value  of  the  angle  9  for  which  the  angular  function  Sm^  (c*cos0) 

is  calculated. 

ARGG:  Temporary  holding  array  for  outputing  ARG. 

ARGI:  Input  parameter  and  initial  value  of  the  angle  9  for  which  the 

angular  function  S^^(c,cos9)  is  desired. 

BLIST:  Array  used  in  the  Bouwkatnp  procedure.  It  contains  the  3“ 

coefficients  defined  in  Eq.  (21). 

C:  c. 


Cl: 

CL: 


CLSPAC: 


COEF: 

CSQ: 

DARG: 

DC: 

DEC: 

DEC2 : 


DNUM: 


Input  parameter  and  initial  value  of  c. 

Eigenvalue  A^.  Initially  equal  to  the  approximation  returned 
by  subroutine  "GETEIG",  then  converges  to  A^  during  Bouwkamp 
procedure. 

The  estimated  spacing  between  two  eigenvalues.  Used  to  determine 
if  the  Bouwkamp  procedure  has  converged  to  the  proper  eigenvalue 
and  to  approximate  an  upper  bound  on  a  new  estimate  if  a  new 
starting  value  is  needed. 

Term  used  in  computing  the  normalizing  factor. 

9 

c~. 

Input  parameter  and  step  size  used  to  generate  ARG. 

Input  parameter  and  step  size  used  to  generate  c. 

A  constant  set  equal  to  10  (^DEC+1).  yse(j  tQ  determine 

convergence  of  the  angular  function  series  of  Eq.  (12). 

-(NDE'C-1) 

A  constant  set  equal  to  10  .  Used  to  determine 

convergence  of  the  Bouwkamp  eigenvalue  procedure. 

The  normalizing  factor  used  to  provide  Meixuer-Schiifke 
normalization  for  the  angular  functions.  It  is  equal  to  the 
normalized  value  of  dn 
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EIG2,  EIG3, 
and  EIG4 : 
ENR: 


EX: 


FL: 

FTERM: 


GLIST: 

IBLIM: 


IM: 


IW6: 

IX: 

IXX: 

JHI: 


L: 

LAM1-LAM5 : 


LIM1 : 
LNUM: 


M: 

MAXAC: 


Previous  eigenvalues  used  to  generate  the  eigenvalue  estimate 
when  £-m  is  large  enough. 

The  array  used  to  contain  the  D  constant  ratios  returned  from 

"CONVER".  The  coefficients  N111  are  calculated  from  these  ratios 

n 

in  the  main  program. 

ENR(I)  -  ^21^21-2  *  is  even» 

ENR(I)  *  d2l+l/d2I-l *  if  iS  °dd* 

Constant  set  to  10^^  .  Used  in  testing  numbers  to  prevent 

overflow  on  the  computer. 

The  eigenvalue  approximation  returned  by  "GETEIG". 

The  largest  term  in  the  series  of  Eq.  (12)  used  to  form  the 
angular  function.  It  is  used  to  estimate  the  resulting 
subtraction  error  present  in  the  angular  function. 

Array  used  in  "CONVER"  as  part  of  the  Bouwkamp  procedure.  It 
contains  the  coefficients  as  described  in  Eq.  (20). 

The  number  of  terms  used  in  the  angular  function  series  of 
Eq.  (12),  determined  as  follows:  IBLIM  »  LIMI/2  -  IX. 

Input  parameter  and  the  increment  used  to  generate  desired 
values  for  m. 

C£-m)/2,  truncated  to  an  integer. 

Equals  zero  if  Z= m  is  even;  equals  one  if  £*m  is  odd. 

IX- 1. 

The  number  of  Legendre  functions  to  be  calculated  for  a  given 
argument,  determined  as  follows:  JHI  *  2*(LNUMfCMAX+NDEC) , 
where  LMAX  =  Cl  +  (NC-1)*DC  is  the  largest  value  desired  for  c. 
1. 

Coefficients  used  in  "GETEIG"  to  generate  the  power  series 
expansion  for  the  eigenvalue  approximation  in  terms  of  c. 
2*(L-M+C+NDEC) . 

The  number  of  successive  values  of  Z  starting  with  Z  *  m  for 
which  angular  function  values  are  desired, 
m. 

The  number  of  decimal  digits  in  the  printed  output  for  the 
angular  functions.  This  parameter  is  set  to  eight  in  this 
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MMIN: 

MNUM: 

NACC: 

NARG: 

NC: 

NDEC: 

HEX: 

P: 

PI: 

PLEGl : 

PNORM: 
PTEMP : 

PTEST: 

R: 


version  of  the  program.  See  the  section  entitled  ACCURACY  OF 
RESULTS  for  information  on  when  and  how  to  change  the  parameter. 
Input  parameter  and  starting  value  for  m. 

Input  parameter  indicating  the  number  of  values  of  ra  for  which 
angular  function  values  are  desired. 

Array  containing  a  measure  of  the  number  of  decimal  digits  that 
are  accurate  in  the  printed  value  for  the  angular  function. 

Input  parameter  indicating  the  number  of  values  of  0  for  which 
angular  function  values  are  desired. 

Input  parameter  indicating  the  number  of  values  of  c  for  which 
angular  function  values  are  desired. 

Initialization  parameter  that  is  set  equal  to  the  number  of 
decimal  digits  available  on  the  user's  machine  in  double 
precision  arithmetic. 

Initialization  parameter  that  is  set  equal  to  the  maximum 
exponent  size  that  is  available  on  the  user's  machine  in  double 
precision  arithmetic. 

A  doubly  dimensioned  array  that  contains  the  ratios  of 
successive  associated  Legendre  functions,  where  P(k,l)  =*  1, 

P(k,j)  *  ^rn+j/^nH-J-l’  witb  Pn  §iven  The  index  k 

refers  to  the  value  of  9.  The  special  cases  of  0  *  0,  90,  and 

180°  are  handled  somewhat  differently,  as  described  in  the 
section  entitled  COMPUTATIONAL  PROCEDURES. 


Value  for  it,  specified  to  36  digits  but  truncated  to  NDEC 
digits  by  the  computer. 

Vector  of  length  NARG  containing  scaling  coefficients  used  to 
prevent  an  overflow  while  forming  P^Cn). 


Equal  to  log^Q  of  P™(n)  as  given  in  Eq.  (16). 

Vector  of  length  NARG  which  contains  values  for  P^(n).  These 
values  are  scaled,  if  necessary,  to  prevent  computer  overflow, 


,m. 


by  10PLEG1. 

Constant  set  to  10  degrees.  Input  values  of  0  are  set  equal 
to  0,  90.  or  180°  if  they  are  within  PTEST  of  these  values. 
■£-m. 
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RL: 

RL2 : 

RM: 

RM2 : 

S: 

SIGN: 

SI: 

DIMENSIONS 


l. 

21. 

m. 

2m. 

A  temporary 
output. 

The  sign  of 
The  angular 

ID  STORAGE 


holding  array  for  the  angular  function  prior  to 


(1) 

function  (c,n). 


The  storage  requirements  for  the  program  PANGFN  are  dominated  by 
dimensioned  arrays.  Everything  else  takes  about  10,000  words  of  storage. 
The  minimum  dimension  requirement  (M.D.R.)  for  each  array  is  determined  by 
the  desired  range  of  parameters  as  follows: 

1.  The  M.D.R.  for  BLIST,  GLIST,  and  ENR  is  given  by 
(LNUM  +  GMAX  +  NDEC)  where  LNUM  is  the  number  of  values 
of  t  desired,  CMAX  *  Cl  +  (NC-1)  *DC  is  the  largest  value 
of  c  desired,  and  NDEC  is  the  number  of  decimal  digits  in 
double  precision  words  on  the  user's  computer.  This 
dimension  is  set  at  250  in  the  listed  version  of  PANGFN. 

2.  The  M.D.R.  for  PLEG1,  PNORM,  and  PTEMP  is  NARG,  the 
number  of  values  of  0  for  which  angular  function  values 
are  desired.  This  dimension  is  set  at  10  in  the  listed 
version  of  PANGFN.  It  can  be  increased  (or  decreased)  if 
more  (or  fewer)  values  of  0  are  desired. 

3.  The  M.D.R.'s  of  the.  doubly-dimensioned  array  P(K,J)  are 
NARG  for  K  and  JHI  =  2* (LNUM  +  CMAX  +  NDEC)  for  J.  The 
value  for  JHI  is  just  twice  that  of  the  M.D.R.  for  BLIST, 

GLIST,  and  ENR  given  above  in  Item  1.  JHI  is  set  equal 

to  500  in  the  listed  version  of  PANGFN. 

4.  All  other  arrays  have  a  dimension  of  three.  This  is 

required  for  the  printed  output  format. 

v. 

If  adequate  computer  storage  is  available,  it  is  advisable  to  set  the 
dimensions  large  enough  to  accommodate  any  anticipated  parameter  input  and 
then  forget  about  them. 


12 


PARAMETER  INPUT 

The  input  to  PANGFN  consists  of  a  series  of  data  cards  as  follows: 

Data  Card  1:  Format  15  -  This  card  contains  the  integer  MMIN,  located  in 

the  first  five  spaces  of  the  card,  right  justified.  MMIN  is 

the  smallest  value  of  m  to  be  used  in  the  computation  of 
the  angular  function. 

Data  Card  2:  Format  15  -  This  card  contains  the  integer  IM,  located  in  the 
first  five  spaces  of  the  card,  right  justified.  IM  is  the 
increment  used  to  generate  subsequent  values  of  m  from  MMIN. 

Data  Card  3:  Format  15  -  This  card  contains  the  integer  MNUM,  located  in 

the  first  five  spaces  of  the  card,  right  justified.  MNUM  is 

the  number  of  values  of  m  for  which  angular  function  values 
are  desired. 

Data  Card  4:  Format  15  -  This  card  contains  the  integer  LNUM,  located  in 

the  first  five  spaces  of  the  card,  right  justified.  LNUM  is 

the  number  of  values  of  •£  for  which  angular  function  values 
are  desired. 

Data  Card  5:  Format  D20.I0  -  This  card  contains  the  value  of  ARGI, 

located  in  the  first  twenty  spaces  of  the  card.  ARGI  is  the 
initial  value  for  9  and  is  used  with  DARG  to  generate  all 
the  desired  values  of  9. 

Data  Card  6:  Format  D20.10  -  This  card  contains  the  value  of  DARG,  located 

in  the  first  twenty  spaces  of  the  card.  DAR  is  the 
increment  used  to  generate  subsequent  values  of  9  from  ARGI. 

Data  Card  7:  Format  15  -  This  card  contains  the  integer  NARG,  located  in 

the  first  five  spaces  of  the  card,  right  justified.  NARG 

is  the  number  of  values  of  9  for  which  angular  function 
values  are  desired. 

Data  Card  8:  Format  D20.I0  -  This  card  contains  the  value  of  Cl,  located 

in  the  first  twenty  spaces  of  the  card.  Cl  is  the  initial 
value  of  c  used  with  DC  to  generate  subsequent  values  of  c. 

Data  Card  9:  Format  D20.10  -  This  card  contains  the  value  of  DC,  located 

in  the  first  twenty  spaces  of  the  card.  DC  is  the  increment 
used  to  generate  subsequent  values  of  c. 
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Data  Card  10:  Format  15  -  This  card  contains  the  integer  NC,  located  in 
the  first  five  spaces  of  the  card,  right  justified.  NC  is 
the  number  of  values  of  c  for  which  angular  functions  values 
are  desired. 

The  program  can  easily  be  modified  if  the  user  desires  to  specify  values 
for  n  *  cos8  rather  than  9.  The  following  changes  are  required: 

1.  Change  ARG  in  statement  72  in  the  main  program  to  BARG. 

Add  BARG  to  double  precision  list. 

2.  Add  statement  ARG  *  DARC0S(ARG*PI/180.DQ)  immediately 
following  statement  72  in  the  main  program. 

3.  Change  statement  150  in  the  main  program  to  read  ARGG(ISTEP)  * 
BARG. 

4.  Change  statment  10  in  subroutine  "PLEG"  to  read  BARG  *  ARG. 

5.  Add  statement  ARG  *  DARCOS(ARG*PI/180.DO)  immediately 
following  statement  10  in  "PLEG". 

6.  Change  the  print  format  for  A(l)  =  n  in  statement  1  of 
subroutine  "OUTPUT"  from  F8.3  to  F8.5  to  provide  five  digits 
to  the  right  of  the  decimal  in  the  printed  value  for  n. 

PARAMETER  RANGES 

PANGFN  was  developed  and  tested  on  the  PDP-11/45  computer  at  the 
Underwater  Sound  Reference  Detachment  fUSRD)  of  NRL  for  the  following 
parameter  ranges: 

0°  *  9  £  180° 

0.00001  <  c  £  100 

0  4  m  4  100 

m  <  l  <  m+100 

PANGFN  is  not  limited  to  these  ranges,  however.  They  were  chosen  to  be 
compatible  with  the  relatively  small  core  memory  of  the  PDP-11/45  computer 
at  the  USRD.  By  increasing  the  dimension  specifications  above  those  given 
in  the  program  listing  in  Appendix  B  to  allow  for  more  terras  in  the  series 
used  to  calculate  the  angular  functions,  the  ranges  on  c  and  £-m  can  be 
increased  indefinitely.  The  minimum  dimension  specifications  required  for 
larger  values  of  c  and  £-m  are  given  in  the  section  entitled  DISTENSIONS  AND 
STORAGE.  There  are  no  limitations  on  the  range  for  m.  Increasing  m  does 
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not  require  any  changes  in  the  dimension  specifications.  Values  of  c  smaller 
than  0.00001  can  be  chosen,  if  desired.  However,  the  output  format 
specification  for  c  given  in  statement  1005  must  be  changed  from  F15.5  to 
include  more  digits  to  the  right  of  the  decimal. 

The  version  of  PANGFN  listed  in  Appendix  B  was  also  run  on  the  Texas 
Instrument  (TI)  ASC  computer  at  NRL  and  tested  for  the  range  of  parameters 
given  above.  The  results  were  consistent  with  those  for  the  PDP-11/45. 

Since  the  TI  ASC  computer  has  more  than  500,000  words  of  core  memory,  the 
dimension  specification  can  be  increased  substantially  on  this  machine.  As 
an  example,  the  dimensions  of  BLIST,  GLIST,  and  ENR  were  increased  to  2500 
and  the  second  dimension  of  P  was  increased  to  5000  (JHI  *  5000) .  Values 
of  the  angular  function  were  then  successfully  computed  for  c  and  -£-m  both 
larger  than  1000. 

COMPUTATIONAL  PROCEDURES 

There  are  four  major  computations  in  PANGFN:  1)  determination  of 
the  eigenvalues,  2)  determination  of  a  normalizing  factor  for  the  D 
constants,  3)  determination  of  the  Legendre  functions,  and  4)  calculation 
of  the  angular  function. 


Determination  of  the  Eigenvalues 

Accurate  eigenvalues  A^  are  obtained  by  use  of  the  variational 
procedure  developed  by  Bouwkamp.  This  procedure,  found  in  subroutine 
"CONVER",  takes  an  approximate  starting  value  for  A^  and  produces  a 
correction  oA^.  This  correction  is  added  to  the  starting  value  to  obtain 
a  better  approximation  to  the  eigenvalue  and  the  Bouwkamp  procedure  is 
repeated  with  this  new  approximation  as  the  starting  value.  Convergence 
of  the  procedure  is  obtained  when  the  relative  contribution  of  the 
correction  becomes  less  than  10  ^  >  where  NDEC  is  the  number  of 

decimal  digits  used  in  the  calculation. 

The  key  to  obtaining  the  correct  eigenvalue  lies  in  the  choice  of  the 


,(D 


The  Bouwkamp  procedure  will  always  converge  to 


initial  approximation  A  ?  . 

me.  n ) 

an  eigenvalue,  but  it  will  only  converge  to  A^  when  A^  is  sufficiently 

close.  Otherwise  it  will  converge  to  another  eigenvalue  A^,  for  the  same 

value  of  m  and  c.  The  reason  for  this  is  that  the  Bouwkamp  procedure  does 


not  depend  explicitly  on  £  but  only  implicitly  through  the  eigenvalue 
A^,  In  addition.,  since  m  is  restricted  to  even  or  odd  values  depending 
on  whether  £-m  is  even  or  odd,  respectively,  Eq.  (17)  has  two  distinct 
forms — one  for  even  £-m  and  the  other  for  odd  £-m.  Therefore,  the  Bouwkamp 
procedure  always  converges  to  a  characteristic  value  such  that  £'  has 

the  same  parity  as  £. 

The  eigenvalues  A^,  £  ■  m,  rn+1,  ...  form  a  monotonically  increasing 
sequence  of  positive  real  numbers.  For  each  eigenvalue  there  exists  a 
region  of  convergence  for  which  the  Bouwkamp  procedure  will  converge 
to  that  eigenvalue.  If  the  value  of  A^  is  slightly  greater  than  the  upper 
bound  of  J2  £  or  slightly  lower  than  the  lower  bound  of  convergence  will 

be  to  V  £+2  °r  Am,£-2’  resPectively- 

Several  different  methods  are  used  to  obtain  approximations  to  the 

eigenvalue  for  starting  the  Bouwkamp  procedure.  These  include  a  power 
series  expansion  in  c2,  an  asymptotic  expansion  in  1/c,  extrapolation  of 
previous  eigenvalues,  and  the  approximation  A  £  *  £(£+1).  The  choice  of 
methods  depends  on  the  parameters  c,  m,  and  £.  Table  I  summarizes  the 


choices. 


Table  I  -  Choice  of  method  to  obtain  eigenvalue  approximation 


£  =*  m,  m+1 


£  *  m+2,  m+3 


£  3  m+4 


£  >  m+4 


c  <  6 

6  <  c  <.  8 

6  <  c  <  8 

c  >  8 

c  >  8 

c  £  5 

5  <  c  i  6 

6  <  c  -  8 

6  <  c  -  8 

c  >  8 

c  >  8 

c  <  8 

c  >  8 

c  >  8 


all 
m  <  4 

m  -  4 

m  <  10  +  c 
m  Z  10  +  c 

all 
all 
ra  <  4 
m  2  4 
m  £  6 
ra  >  6 

all 
m  <  3 
m  2-  3 


Method 

c2  expansion 
I/c  expansion 
c“  expansion 
I/c  expansion 
£(£+1)' 

c2  expansion 
extrapolation 
1/c  expansion 
extrapolation 
i/c  expansion 
extrapolation 

extrapolation 
1/c  expansion 
extrapolation 

extrapolation 
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This  method  does  not,  however,  guarantee  that  lies  in  3^;  it 

simply  guarantees  that  is  somewhat  close  to  A^.  For  some  m,  •£,  and  c 

the  region  A^  may  be  extremely  narrow  and,  therefore,  the  estimate  A^^ 

would  have  to  be  very  close  to  A  n  to  achieve  convergence.  No  simple 

m£  (1) 

method  can  absolutely  guarantee  that  A  £  lies  in  £2^.  A  procedure  has 
therefore  been  developed  to  determine  if  the  resulting  value  A^  is  actually 
the  correct  eigenvalue.  In  order  for  A  £  to  be  consistent  with  the  other 
eigenvalues  previously  obtained,  it  must  satisfy  the  following  conditions: 


>  ^ra,-£-l’ 


A  }  -  A(P  <  A(P  -  A  f  • 
m-c  m-c  me  m,-e-l 

If  A  n  passes  both  of  these  tests,  it  is  the  correct  eigenvalue.  If  A^ 
fails  the  first  inequality,  it  is  equal  to  A  n  ,  In  this  case 

/  n  ra,e.-4. 

A  can  be  considered  as  a  lower  bound  for  A  ».  An  upper  bound  can  be 
me  qv  m-e 

determined  by  taking  A^?  +  1.5*CLSPAC  where  CLSPAC  is  the  previous  eigenvalue 


spacing  (Am>£_1  -  Am>£_2). 


If  A^  fails  the  second  inequality,  it  is  equal  to  A^  £+2« 


In  this 


case  can  be  considered  as  an  upper  bound  for  A^.  A  lower  bound  can  be 

taken  as  A  *>  .  .  Once  upper  and  lower  bounds  for  the  eigenvalue  have  been 
m,-u- 1 

established,  a  new  starting  value  is  taken  as  the  mean  of  the  bounds.  The 

Bouwkamp  procedure  is  then  repeated  with  this  new  approximation.  The 

resulting  eigenvalue  is  tested  to  see  whether  it  lies  within  the  hounds. 

If  so,  has  been  obtained  and  the  process  is  concluded.  If  the  resulting 

eigenvalue  is  less  than  the  lower  bound,  then  A^  must  lie  between  the 

starting  value  and  the  upper  bound,  so  that  the  starting  value  becomes  the 

new  lower  bound.  If  the  resulting  eigenvalue  is  greater  than  the  upper 

bound,  then  A  0  must  lie  between  the  lower  bound  and  the  starting  value, 
mc 

so  that  the  starting  value  becomes  the  new  upper  bound.  A  new  starting 
value  is  taken  to  be  the  mean  of  the  new  bounds,  and  the  process  is 
repeated  until  convergence  to  A^£  is  obtained. 
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Determination  of  the  Normalizing  Factor 

The  variational  procedure  of  Bouwkamp  includes  as  an  intermediate  step 

the  calculation  of  the  ratios  *  ct  d  /d  „.  Dividing  by  cx  as  given  by 

Eq.  (22)  results  in  the  values  of  d  /d  These  ratios  are  well-behaved 

ix  n—z 

throughout  the  required  range  of  n.  They  do  not  become  either  extremely 

large  or  extremely  small.  It  is  convenient,  for  computational  purposes, 

to  set  the  largest  D  constant  equal  to  unity.  The  largest  value  occurs 

in  the  region  of  n  *  £-m,  therefore,  d»  is  set  equal  to  one.  In  theory 

all  of  the  other  D  constants  can  be  obtained  from  d0  bv  successive 

C-m  ' 

multiplication  of  the  ratios  d^/d^.  In  practice,  however,  these 
coefficients  are  not  evaluated  directly  in  PANGFN  due  to  the  computer 
underflow  that  would  likely  result  from  their  extreme  range  in  magnitude. 

Angular  function  values  obtained  from  Eq.  (12)  with  dy  set  equal 
to  unity  require  a  normalizing  factor  to  provide  the  desired  Meixner- 
Schafke  normalization.  The  first  step  in  calculating  this  factor  is  to 
evaluate  the  left-hand  side  of  Eq.  (24) ,  starting  at  n  *  £-m  with  d^  ^  set 
equal  to  unity  and  proceeding  both  upward  and  downward  in  n  until  convergence 
results.  Each  term  in  the  series  is  obtained  from  the  previous  term  by 
multiplication  (upward)  or  division  (downward)  by  both  the  square  of  the 
appropriate  D  constant  ratio  and  by  a  ratio  of  integers  to  account  for  the 
coefficients.  Convergence  is  assumed  when  the  relative  contribution  of  the 
last  term  (both  upward  and  downward)  is  less  than  10  (NWJC+1)^  where  ^DEC  is 
the  number  of  decimal  digits  in  double  precision  words  in  the  computer. 

This  procedure  should  prevent  underflow  or  overflow  during  evaluation  of  the 
series.  The  normalising  factor  DNUM  is  then  obtained  by  dividing  the 
resulting  sura  into  the  right-hand  side  of  Eq.  (24)  and  taking  the  square  root. 
This  quantity  is  the  magnitude  of  the  normalized  value  for  d£_m  required 
to  provide  Meixner-Schafke  normalization  of  the  angular  functions.  The 


algebraic  sign  of  dp  is  obtained  by  progressive  multiplication  of  the 

-c-ra 

signs  of  the  ratios  of  the  D  constants,  beginning  with  d^/d^  or  d^/d^ 

depending  on  whether  £-m  is  even  or  odd,  and  continuing  to  d^^/d^ 

Since  dg  and  d^  are  always  positive,  this  product  of  signs  is  equal  to  the 

sign  of  d,_m#  The  angular  functions  are  then  obtained  by  multiplying  the 

sum  of  Eq.  (12)  by  dp 

c-m 
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Determination  of  the  Associated  Legendre  Functions 

The  associated  Legendre  functions  required  in  Eq.  (12)  are  calculated 
in  the  subroutine  "PLEG".  Actually,  ratios  of  successive  associated 
Legendre  functions  are  calculated  instead  of  the  functions  themselves  because 
of  possible  overflow  problems  when  m  is  large.  The  recursion  relation  of 
Eq.  (14)  is  modified  for  this  purpose  to  give: 

T%)  =  P®(n)/P  ,(n)  -  ,  n  >  nrf-1  ,  (25) 

n  n  n-l  n-m  (n_m)T* 

n-i 

with  T111  *  1  and  T  ™  *  (2m+l)n.  The  ratio  Tm  is  well  behaved  for  n  j  0; 
ra  nrrl  n 

i.e.,  for  9  ^  90°.  If  0  *  90°,  however,  is  equal  to  zero  for  n-m  odd, 
and  T^  is  unbounded  or  equal  to  zero,  depending  on  whether  n-m  is  even  or 
odd.  The  relevant  ratios  in  this  case  are  those  of  successive  nonzero  values 
of  the  associated  Legendre  functions.  These  ratios  are  calculated  using 

T^O)  *  Pm(0)/P  m„(0)  *  -(n+m-l)/n-m  ,  n-m  even  ,  (26) 

n  n  n-2 

with  ^(O)  -  1. 

m  m 

The  ratios  T  (0)  for  n-m  even  are  stored  in  the  odd  J  locations  of 
n 

P(K,J).  For  convenience  in  calculating  the  successive  terms  in  Eq.  (12), 

the  ratios  1^(0)  f°r  n-m  odd  are  set  equal  to  unity  and  stored  in  the  even 

J  locations  of  P(K,J).  To  avoid  possible  numerical  difficulties  occurring 

near  9  *  90°,  the  program  sets  9  equal  to  90®  if  it  is  within  PTEST  «  10  ' 

degrees  of  this  value.  Th<?  constant  PTEST  is  defined  near  the  beginning 

of  the  main  program.  The  value  of  PTEST  can  be  decreased  if  the  user 

desires  angular  function  values  for  9  closer  than  10  '  degrees  to  90®. 

The  associated  Legendre  functions  and  thus  the  angular  functions  are 

equal  to  zero  when  n  “  ±1;  i.e.,  when  0  *  0  or  180®  and  m  is  unequal  to 

—7 

zero.  If  the  input  value  for  9  is  within  PTEST  *  10  degrees  of  0  or  180° 
and  if  ra  is  unequal  to  zero,  the  program  assumes  that  the  argument  is  equal 

to  0  or  180®,  respectively,  bypasses  the  calculation  of  'Im,  and  directly 

(1)  m 

sets  the  angular  function  S  £  (c,±l)  equal  to  zero. 


For  purposes  of  calculating  the  angular  functions  S^£^(c,n)  using  the 

series  of  Eq.  (12),  it  is  convenient  to  sec  P£(n) ,  the  associated  Legendre 

function  that  is  multiplied  by  the  D  constant  d^  ,  equal  to  unity,  The 

largest  term  in  the  series  is  now  equal  to  or  near  unity.  This  choice 

eliminates  potential  underflow  and  overflow  problems  in  the  series.  The 

sum  of  the  series  must  then  be  multiplied  by  P£(n)  to  correct  for  this  choice. 

The  required  value  for  P^(n)  is  obtained  from  P®  ^(n),  which  was  used  in 

the  same  way  for  ,(c,n),  by  multiplication  by  the  stored  ratio  P;>(n)/ 

P1?  ,(n).  The  value  for  Pra(n) ,  required  for  ^(c,n),  is  calculated  using 
■L—  i  m  mm 

Eq.  (16).  Since  PjSkn)  can  be  extremely  large  when  m  is  large,  the  value  for 
P^(n)  is  calculated  and  stored  as  a  logarithm  to  the  base  10  to  avoid 
overflow. 

Determination  of  the  Angular  Function 

The  angular  functions  are  calculated  using  Eq.  (12).  The  series  is 
evaluated  by  starting  at  n  *  £-m  and  proceeding  both  upward  and  downward 
in  n  until  convergence  is  obtained.  Each  term  in  the  series  is  derived  from 
the  preceding  te:.m  by  use  of  successive  multiplication  (upward)  or 
division  (downward)  of  the  ratios  of  D  constants  and  the  associated 
Legendre  functions.  To  start  the  process,  both  d^  ^  and  P^n)  and  thus  the 
initial  term  dff  P^(n)  are  set  equal  to  unity.  The  logarithm  to  the  base 
ten  of  the  resulting  sura  is  now  taken  and  added  to  the  corresponding 
logarithm  of  P°(n)  and  d^_m.  This  corrects  for  setting  both  P^(n)  and 
d^  ^  equal  to  unity.  The  resulting  logarithm  of  the  angular  function  is 
passed  to  the  subroutine  "OUTPUT".  "OUTPUT"  separates  the  mantissa  from 
the  characteristic,  takes  its  antilogarithm,  and  combines  it  with  the  proper 
sign  to  obtain  the  base  of  the  angular  function.  The  number  is  the  output 
in  two  parts:  the  base  (-9.9999999  £  base  <  9.9999999)  and  its  exponent 
(-999  £  exponent  £  999)  as  given  by  the  characteristic.  If  the  angular 
function  becomes  as  large  as  101000.  possibly  for  very  large  m,  or  as  small 
as  10  ,  possibly  for  9  £  10  degrees  or  (180-9)  £  10  degrees,  then  the 

exponent  printout  must  be  appropriately  expanded. 
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COMPUTATION  TIME 


The  execution  for  PANGFN  time  depends  on  the  input  data.  In  particular, 
larger  values  of  c  and  £-m  take  more  time.  The  following  examples  are 
representative  of  the  execution  times  for  PANGFN  on  the  PDP-11/45  computer 
at  USRD  and  the  TI  ASC  computer  at  NRL  for  9  *  0®,  10°,  90° t  m  «  0; 

Z  *  0,  1,  100;  and  selected  values  of  c  (see  Table  II).  The  execution 

times  are  nearly  independent  of  m;  e.g.,  the  same  times  as  shown  here  are 
obtained  for  ra  *  100  and  Z  *  100,  101,  ...,  200.  The  execution  times  where 
only  one  value  of  9  is  desired  are  greater  than  half  of  those  shuon  above. 
Therefore,  it  is  economical  to  include  all  desired  values  of  9  in  a  single 
run.  (If  more  than  10  values  of  9  are  desired,  some  of  the  dimension 
specifications  must  be  increased.  See  the  section  entitled  DIMENSIONS  AND 
STORAGE.) 


Table  II  -  Execution  time  for  selected  values  of  c 


c  PDP-11/45  TIME  TI  ASC  TIME 


1.0 

42  s 

2.0  s 

10.0 

48  s 

2.4  s 

50. C 

70  s 

3.2  s 

100.0 

96  s 

4.3  s 

PRINTED  OUTPUT 

The  output  from  PANGFN  consists  of  numerical  tables,  as  shown  in 

Appendix  A.  Numerical  values  for  the  eigenvalue  A  py  the  desired  arguments 

m-c  ( i ) 

9  together  with  the  corresponding  angular  functions  (c,cos6),  and 
accuracy  estimates  are  given  for  desired  choices  of  c,  m,  and  Z. 

The  argument  (ARG)  is  printed  with  three  digits  to  the  right  of  the 
decimal  point.  The  eigenvalues  (EIG)  and  the  angular  functions  (S)  are  each 
printed  with  eight  decimal  digits.  The  accuracy  estimate  (ACC)  is  printed 
as  an  integer  and  indicates  the  number  of  leading  decimal  digits  in  the 
printed  output  that  are  likely  to  be  accurate.  A  discussion  of  the  accuracy 
estimate  is  given  in  the  section  entitled  ACCURACY  OF  RESULTS.  A  procedure 
for  changing  the  number  of  printed  decimal  digits  in  the  output  is  outlined 
under  the  following  section. 
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ACCURACY  OF  RESULTS 


The  procedure  used  to  obtain  the  prolate  eigenvalues  is  well-behaved 
numerically.  Therefore  all  eight  digits  printed  in  the  output  will  be 
accurate,  with  the  last  digit  rounded. 

The  series  expansion  used  to  evaluate  the  angular  function  is  not  always 
well-behaved  numerically.  Subtraction  errors  can  occur  in  the  summation  of 
the  series,  especially  for  large  c  and  low  m  and  £.  A  good  approximation  for 
the  number  of  decimal  digits  of  accuracy  lost  to  subtraction  error  in 
this  summation  is  given  by  SE  *  LOG  ( i FTERM/S 1 [ )  where  FTERM  is  the  largest 
term  in  the  series  and  Si  is  the  sum  of  the  series.  It  is  estimated  that 
roundoff  error,  slight  inaccuracies  in  the  D  constants,  and  the 
representation  of  the  angular  function  as  a  logarithm  may  contribute  an 
additional  loss  of  two  decimal  digits  of  accuracy.  Since  NDEC  is  the  number 
of  decimal  digits  available  for  double  precision  words,  T  =  NDEC-SE-2  is  the 
expected  accuracy  of  the  calculated  value  of  the  angular  function.  If  T 
is  greater  than  eight,  it  is  reduced  to  eight  to  correspond  to  the  number 
of  decimal  digits  printed  in  the  output.  If  T  is  less  than  zero,  it  is 
set  equal  to  zero.  The  resulting  value  for  T  is  stored  in  the  array  NACC(K) 
and  output  under  the  heading  ACC. 

The  only  time  that  an  accuracy  less  than  eight  decimal  digits  is  likely 
to  be  encountered  is  when  c  is  larger  than  about  30,  and  when  m  and  £  are 
somewhat  less  than  c.  The  subtraction  error  of  SE  digits  obtained  in  this 
case  corresponds  to  a  value  of  S^^(c,n)  near  10  ^  P^(n),  the  largest 

term  in  the  series  used  to  generate  5^£^(c,n)  having  a  value  near  P^(n) .  In 


other  regions  of  the  parameters  c,  m,  and  £,  the  value  of  the  angular 
function  is  usually  accurate  to  all  eight  digits. 

If  the  user  wishes  to  change  the  number  of  decimal  digits  printed  in 
the  output  of  this  program  from  eight,  he  should  make  the  following  changes: 
L.  Change  the  Format  statement  numbered  1  in  subroutine 
"OUTPUT"  to  the  desired  number  of  digits. 

2.  Change  the  Format  statements  numbered  1009,  1010,  and 
1011  in  the  main  program  to  line  up  the  headings  to 
correspond  to  change  1  above. 

3.  Change  the  value  of  MAXAC  to  the  desired  number  of 
printed  digits. 
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4.  Change  statement  number  10,  the  statement  four  lines 
before  statement  number  10,  and  statement  number  200  in 
the  subroutine  "OUTPUT"  to  contain  MAXAC  nines  (9's); 
i.e.,  as  many  as  there  are  decimal  digits  in  the  output. 
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appendix  a 

OUTPUT  EXAMPLE 


10.00000 


S  EIG=  0.355880860  02 

ft  KG  g 

0,000  0 .00000000+000 

30.000  3 • 39  484610+000 

£0,000  3 .63758710+002 
90,000  1*11251 180+003 


ARG 

0.000 

30.000 

40.000 

90.000 


EIG=  0.574504850  02 

S 

00  0 . 00000000+000 
00  5.02814740+001 
00  1.98247900+003 

00  0,00000000+000 


0,796°3327D  02 

0,000  o.oooooooo+ooo  A8C 

^°,00°  3*71580550+002  8 

90  one  e  '^^135040  +  003  | 

0.000  -U. 44030440+003  8 

®  EIG=  0  *  1017121 10  03 

ftKG  g 

0,0000°000+000  A8C 

30.000  1.81584280+003  g 

|0, 000  1*10081440  +  004  8 

90.000  0.00000000+000  g 


9  EIG: 
ARG 
0.000 
30.000 
40.000 
90.000 


!  0.124293780  03 

S 

0.00000000+000 
6 . 58232720+003 
4 . S407442D+003 
1 ■ 64884440+004 


ARG 

10.000 

40.000 

70.000 


7*44511370-003 
2.07418930+001 
5 . 85204700+002 


ARG 

20.000 
50 . 000 
80.000 


9 

3.14874790-001 
8*49424410+001 
9 . 44984710+002 


10  £IG 
ARG 
0.000 
30.000 
40.000  • 
90.000 


=  0*147478230  03 

S 

0.00000000+000 
1  *  88557140+004 
■2.04382750+004 

o.oooooooo+ooo 


ARG 

10.000 

40.000 

70.000 


t  ACC  c 

2.59053280+00*  g  so‘000  5  *  29323350+000  ^ 

•*  8:SS  !:2£{SR{«  i 


ARG 

10.000 

40.000 

70.000 


3 

1.27318050+000 

1.59202200+003 

4.85143980+003 


ARG 

i0,000  7*7742883D+000 
7°,00°  4 . 33944870+003 
70.000  -1.69038410+003 


ARG  s 

I0,000  3.57271300+001 
70,000  1*81417870+004 
70.000  -1.72677970+004 


ARG  g 

10.000  1.32418380+002 
70,000  3.91865480+004 
70,000  -1.74178390+004 


ACC  ARG  s 

8  9  ‘  ^^-80310+001 
8  r n‘00°  4 . 16385470+003 
8  80.000  -1.52944230+003 


ARG  e 

- . 50715800+002  T 
jO.OOO  1.22168960+004  ff 

80.000  -1.00928400+004  8 


ARG  e 

|°'000  *.05631020+003  T 
I0,000  2.33032810+004  g 
80.000  -2.82489820+003  f 


ARG  g 

;!•??“  3-57130..D.OO3  T 
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APPENDIX  B 
PANGFN  LISTING 


PROGRAM  PANGFN 
C 

C  A  FORTRAN  COMPUTER  PROGRAM  FOR  CALCULATING  NUMERICAL 

C  VALUES  OF  THE  PROLATE  ANGULAR  FUNCTION  OF  THE  FIRST  KIND 
C  AND  ITS  ASSOCIATED  EIGENVALUES. 

C 

DIMENSION  8LIST(250)»GLIST(250)jENR(250>,P(10»500)»S(3), 

1 ARGG ( 3 ) , I IS  <  3 ) r PNORM ( 10 ) , NACC ( 3 ) ,PLEG1<10) ? PTEMP(IO) 

DOUBLE  PRECISION  AJ t ARG , ARGG , ARG1 > ARR 
DOUBLE  PRECISION  BLIST 
DOUBLE  PRECISION  CLjCOEF 

DOUBLE  PRECISION  DARG  , DEC  ,  DC  , DNUM  , DNEU  » D0LD 
DOUBLE  PRECISION  EA t EIG2 , EIG3 > EIG4 , EIG5 » ENR > EX , EX1 
DOUBLE  PRECISION  FTERM 
DOUBLE  PRECISION  GLIST 
DOUBLE  PRECISION  C,CSQ,C1 

DOUBLE  PRECISION  P , PI , PLEG1 , PNORM , PTEMP » PTEST 
DOUBLE  PRECISION  RM  , RM2  >  RNDEC 
DOUBLE  PRECISION  S,S1»SIGN 
COMMON  /BLK2/  BLIST » GL I  ST » ENR 
COMMON  /BLK3/  DEC  »  NDEC 
C 

C  THE  FOLLOWING  TWO  STATEMENTS  MUST  BE  MODIFIED  TO  CONFORM  WITH  THE 
C  USERS  COMPUTER 

C  NDEC  J  THE  MAXIMUM  NUMBER  OF  DECIMAL  DIGITS  AVAILABLE  ON  THE 
C  USER'S  COMPUTER  IN  DOUBLE  PRECISION  ARITHMETIC. 

C  NEX:  THE  MAXIMUM  EXPONENT  AVAILABLE  FOR  A  DOUBLE  PRECISION 

C  NUMBER 

C 

N  D  E  C  =  1 6 
NEX=75 
C 
C 

P 1  =  3. 14159265350979323846264338327950288 DO 

MAXAC=8 

RNDEC=NDEC+1 

DEC=10.D0**(-RNDEC) 

EX=10.D0**(NEX~5) 

EX  1  =  10 . DOS* ( 5-NEX ) 

PTEST  =  1 . D-7 
C 


READ  IN 

INITIAL  ARGUMENTS 

READ 

1101 , MM IN 

READ 

1101,  IM 

READ 

1101  ,MNUM 

READ 

t  lOlrLNUM 

READ 

1102, ARG1 

READ 

1102, DARG 

READ 

1 101 ,NARG 

* 
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READ  1102, Cl 
READ  1102, DC 
READ  1101, NC 
C 

C  BEGIN  PROGRAM 

C  OUTER  LOOP  M  LOOP 

NEXT  LOOP  C  LOOP 

INSIDE  LOOP  L  LOOP 

DO  500  IDUMM=1 , MNUM 
M=MMIN+(IDUMM-1)*IM 
RM=M 

RM2=2.D0*RM 

JHI=2*(LNUM+Cl  +  (NC-l)*DC-i-NDEC) 

CALL  F'LEG  ( M  ,  JHI ,  ARG 1 ,  DARG  ,  N ARG  ,  P ,  PNORM ,  PI ,  PTEST ) 
DO  400  JNC=1 , NC 
DO  3  K=1 ,NARG 

PTEMP  <  K )  =  1 . DO 
PLEG 1 ( K ) =PNORM  <  K  > 

CONTINUE 
C=Cl+< JNC-1)*DC 
CSQ=C*C 
EIG2=0.D0 
EIG3=0 • DO 
EIG4=0 . DO 
EIG5=0  *  DO 
PRINT  1005, C,M 
DO  300  IL=1 , LNUN 
L=M+IL~1 

GET  THE  STARTING  EIGENVALUE 

CALL  GETEIG(L,M,C,CL,EIG2,EIG3,EIG4,EIG5) 

REFINE  THE  STARTING  VALUE 

CALL  CONVER ( L , M , C , CL , EIG3 , EIG5 > 

IU6=(L-M)/2 
IX=L-M-2*IU6 
IXX=IX-1 

IF(L.EQ.M)  GO  TO  7 
DO  5  K=1 , NARG 

PTEMP  <  K  )  =PTEMP  ( l\  )  *p  (  K ,  L-M  +  l ) 

IF <DABS( PTEMP (K) ) .  LT .EX)  GO  TO  5 
PLEG1 ( K  >  =  DLQG10 ( DABS (PTEMP ( K ) ) HPLEG1 (K) 

IF ( PTEMP (K) . LT ♦ 0 . DO ) PTEMP ( K  >  =- 1  *  DO 
IF ( PTEMP (K) .GT. 0 . DO ) PTEMP ( K ) =1 . DO 
CONTINUE 
CONTINUE 

LIM1=2*(L-M+C4-NDEC> 
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20 

30 

C 

C 

C 


50 

60 


70 

71 


C 

C 

C 

,<c 


IBLIM=LIMl/2-IX 
SIGN=1 .DO 
DO  30  1=1 > IBLIM 
ARR=IX+I+I 
£A=ARR+ARR+RM2 
IF ( I . GT . IU6 )  GO  TO  20 
SIGN=SIGN*(ENR(I)/DABS<ENR(I) ) > 

ENR ( I )  =  ( EA-1 . DO )*  <  EA+1 . BO )*ENR ( I ) / < ( ARR+RM2 ) * 
1  (ARR+RM2-1  *D0)icCSQ) 

CONTINUE 

COMPUTE  THE  NORMALIZING  FACTOR 


1 

1 


1 


DNUM=1 .DO 
C0EF=1 ♦ DO 
JL0U=L-M+2 
JTERM=IU6 

DO  50  J= JLOU »  LIM1 »  2 


A  J  =  J 

JTERM=JTERM+1 

COEF=COEF¥<  <AJtRM2)X<ENRC  JTERM)/AJ)*<  (AJ+RM2-1  .DO)* 
ENR( JTERM)/(AJ-1 .DO) ) * C A J*2 , D0+RM2-3 . DO ) / ( A J*2 . D0+RM2 
■f  1 .  DO  ) 

DNUM=DNUM+COEF 

IF ( DABS ( COEF/DNUM ). LT . DEC )  60  TO  60 
CONTINUE 
JLOU=L-M 

IF( JL0U.LT.2)  GO  TO  71 
C0EF=1 .DO 
JTERM=IU6 
J= JLOU 

DO  70  J J  =  2  » JLOU  >  2 
A  J  =  J 


C0EF=C0EF*<AJ/<AJ+RM2)/ENR( JTERM) )*< ( A J- 1 . DO ) / ( AJ+RM2 


- 1 . DO ) /ENR  <  JTERM ) )*(AJ*2 
JTERM= JTERM- 1 


D0tRM2t 1 . DO ) / ( A  J*2 


DO- 


J  =  J  -  2 

DNUM=DNUM+C0EF 

IF ( DABS ( COEF/DNUM ) . LT . DEC )  GO  TO  71 
CONTINUE 

DNUM=  1  ,  DO/DSQRT  ■:  DNUM  ) 

ISTEP=1 

JL0U=IU6+1 

PRINT  1006  >  L  >  CL 

IF ( NARG » GE . 3) PRINT  100? 

IF ( NARG . EG . 2 ) PRINT  1010 
IF ( NARG . LE . 1 ) PRINT  1011 


COMPUTE  THE  ANGULAR  FUNCTION  S 


DO) 
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DO  200  K=1 t NARG 

ARG=ARGH-(K-1  ♦DO)*DARG 

FOR  ARG=0 » 90 » OR  180  DEGREES  USE  SPECIAL  HETHODS  TO  DETERMINE 
ANGULAR  FUNCTION 

IF ( ( DABS ( ARG-90 . DO ) .LT . PTEST ) . AND . ( IX . EQ . 1 ) )  GO  TO  7 
IF ( (DABS< ARG-180 .DO ) .LT .PTEST) . AND* (M ♦ NE . 0 ) )  GO  TO  7 
IF < (DABS(ARG) .LT. PTEST) .AND. (H.NE.O)  )  GO  TO  75 
GO  TO  80 
75  S ( ISTEP ) =-EX 

I IS ( ISTEP ) =1 
NACC( I STEP )=M AX AC 
GO  TO  150 
80  FTERM=EX1 

D0LD=1 .DO 

IF  (PTEMF'(K)  ♦ LT . 0 , DO )  D0LD=-1 .  DO 
S1=D0LD 

DO  100  J=JLQU » IBLIM 

DNEW=DOLD*ENR(  J )  :KP  (  K  »  J+ J+IXX+2) *P < K »  J  +  J  +  IXX+1) 
S1=S1+DNEU 

IF ( DABS ( DNEU ) . GT . FTERM ) FTERM=DABS ( DNEU ) 

IF ( SI  *  EQ . 0 ♦ DO )  GO  TO  95 
IF ( DABS ( DNEU/S1 ) . LT . DEC )  GO  TO  101 
95  DOLD=DNEW 

100  CONTINUE 

101  IF ( IW6 . LT ♦ 1 )  GO  TO  111 
D0LD=1 ♦ DO 

IF (PTEMP(K) . LT ♦ 0 . DO ) D0LD=-1 . DO 
J  =  IU6 

DO  110  J J= i > IU6 

DNEU=DQLD/(P(K,  J-f  J  + IXX  +  2 > *P ( K >  J  +  J  +  IXX  +  1 )  >/ENR( J) 
S1--S1  +  DNEU 

IF (DABS (DNEU) . GT . FTERM ) FTERM=DABS ( DNEU > 

IF ( SI . EQ . 0 . DO )  GO  TO  109 
IF(DABS(DNEW/S1 ) .LT.DEC)  GO  TO  111 

109  DQLD=DNEU 

J  =  J-1 

110  CONTINUE 

111  IF ( SI . EQ ♦ 0  *  DO  )  GO  TO  120 

S (ISTEP ) =DLOG 1 0 ( DABS ( S 1 ) ) iDLOG 1 0 ( DABS ( DNUM ) ) + 

1  DL0G10 ( DABS ( PTEMP ( K ) > ) FPLEGl (K) 

1 1 S ( ISTEP  >  =  1 

IF (SI .LT.O.DO) IIS( ISTEP >=-l 
IF(SIGN.LT.O.DO) II3( I  STEP ) =- 1 *1 1 S ( I  STEP ) 

GO  TO  125 

120  3 ( ISTEP ) =-EX 

1 1 S ( I  STEP )  =  1 

125  NACC( ISTEP )=DLOG 10 (DABS ( ( FTERM+EX 1 ) / ( S1FEX l ) ) ) 

NACC( ISTEP )=NDEC-N ACC ( ISTEP) -2 

* 
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IF  ( NACC  ( ISTEP  )  »LT  .0)NACC(ISTEP)=0 
IF  ( NACC  ( ISTEF' )  .  GT  .  MAXAC )  NACC  ( ISTEP  )  =MAXAC 
,150  ARGG(ISTEP)=ARG 

C 

C  TIME  TO  OUTPUT 
C 

IF ( ISTEP ♦ NE . 3 )  GO  TO  160 
NUM=3 

CALL  OUTPUT  <ARGG» I I S > S » NACC » NUM ) 

ISTEP=0 

160  ISTEP=ISTEP+1 

IF (K ♦ NE ♦ NARG ♦ OR . ISTEP  «E0  *  1 )  GO  TO  200 
NUM=ISTEP-1 

CALL  OUTPUT<ARGG»IIS»S>NACCiNUM) 

200  CONTINUE 

PRINT  1007 
300  CONTINUE 

400  CONTINUE 
500  CONTINUE 

1005  F0RMAT(43X»2HC=.F15.5>?X>2HM=,I6//) 

1006  F0RMAT(1X»2HL=»I6»6H  EIG=»D15.8> 

1007  FORMAT ( / ) 

1003  F0RMAT(1H1»49X> 23H PRO LATE  ANGLE  FUNCTIONS/) 

100?  F0RMAT<3(8X>3HARG>10X»lHSf9X»3HACCr3X)) 

1010  FORMAT ( 2  <  8X » 3HARG 1 10X  » 1HS  >  9X » 3HACC ,  3X )  ) 

1011  FORMAT  ( 8X  »3HARG>10Xf  lH.S»9Xr3HACC) 

1101  FORMAT ( 15 ) 

1102  FORMAT ( D20 ♦  10 ) 

9999  CONTINUE 

END 
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SUBROUTINE  PLEG  ( H >  JHI »  ARG1  >  BARG  .  NAR6  >  P  >  PNORM » F' I ,  PTEST ) 

DOUBLE  PRECISION  API t ARG>  ARG1 
DOUBLE  PRECISION  BARG 
DOUBLE  PRECISION  DARG 
DOUBLE  PRECISION  P > PI  *  PNORM  ? PTEST 
DOUBLE  PRECISION  R  J  >  RJ2  t RM  >  RMP 
DIHENSICN  P ( 10 » 500 ) » PNORM < 10) 

INITIALIZE  VARIOUS  COEFFICIENTS 

API=PI/180 • DO 
RM=M 

MP=2*M-1 
DO  200  K=1 >NARG 

ARG=AR61+ ( K-l ♦ DO ) *DARG 
BARG=DCOS(ARG*API) 

PNORM (K)=0, DO 

FOR  ARG=0 / 90 i  OR  ISO  DEGREES  USE  SPECIAL  METHODS  TO  DETERMINE 
LEGENDRE  FUNCTION  RATIOS 

IF(DABS<ARG-?0. DO) . LT . PTEST)  GO  TO  150 

IF< <DABS(ARG-130. DO) .LT. PTEST) .AND. (M.NE.O) )  GO  TO  190 
IF< (DABS(ARG) .LT, PTEST) .AND. (M.NE.O) )  GO  TO  190 

NORMAL  COMPUTATION  OF  LEGENDRE  FUNCTION  RATIOS 

0  P ( K » 1 )  =  1 • DO 

P(K>2)=(2,D0*RM+1.D0)*BARG 
DO  100  J  =  3  r  J H I 
R J= J+RM 
RJ2=2.D0*RJ 

P  <  N  j J )  =  ( < R J2-3 . DO ) *BARG- <  R J  +  RM-2 . DO ) /P  <  K »  J- 1 ) ) / < R J-RM- 1 . DO ) 
00  CONTINUE 

IF(M.EQ.O)  GO  TO  200 
I  M=MP 

DO  120  1 1 M= 1 »  MP  »  2 
R  M  P  =  I M 

PNORM (K ) =PNORM ( K) +DL0G1Q ( RMP ) +DLQG10 ( DABS ( DSIN ( ARGfcAPI ) ) ) 
IM=IM-2 

20  CONTINUE 

GO  TO  200 

COMPUTATION  OF  LEGENDRE  FUNCTION  RATIOS  FOR  ARG=90  DEGREES 


150  P ( K  f 1 )  =  1 . DO 

DO  160  J=3,JHIr2 
F:J  =  J+RM 

P  <  K  t J ) =- ( R J+RM-2 . DO ) / <  R J-RM- 1 . DO ) 
P ( K  f  J- 1 )  =  l . DO 
t 
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160  CONTINUE 

IF ( M  *  EQ  » 0  >  GO  TO  200 
IM=MP 

DO  170  1 1 M=1 »  MP »  2 
RMP=IM 

PNORM<K)=PNORh(K)+DLQG10<RMP) 

IM=IM-2 

170  CONTINUE 

GO  TO  200 

190  DO  195  J= 1 f  JHI 

195  P(K » J) =0 ♦ DO 

200  CONTINUE 
RETURN 
END 

* 
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SUBROUTINE  GETEIG < L  , M ,  C  >  CL ,  El 02 » EIG3  ,  EIG4  .•  EIG5 ) 

DOUBLE  PRECISION  BOSH 

DOUBLE  PRECISION  CL, C, CSG 

DOUBLE  PRECISION  EIG2 , EIG3 , EIG4 , EIG5 

DOUBLE  PRECISION  R,RL»RL2,RM,RM2 

DOUBLE  PRECISION  LAM1 , LAM2  , LAM3  > LAM4 

COMPUTE  SOME  HEAVILY  USED  VARIABLES 

CSG=C*C 

RL=L 

RL2=2 • D0*RL 
RM=M 

RM2=2.D0*RM 

R=RL-RM 


THE  CASE  OF  SMALL  M  AND  LARGE  C 

IF ( L ♦ EG ♦ ( M+4 ) *  AND . C ♦ GT . 8 . DO . AND . M . LT .3)  GO  TO  200 
APPROXIMATIONS  BASED  ON  PREVIOUS  EIGENVALUES,  IF  AVAILABLE 
IF ( L . GT . ( M+3 ) )  GO  TO  30 

USE  EXPANSION  IN  TERMS  OF  CSG  FOR  LOU  C,  AND  THE  EXPANSION 
IN  TERMS  OF  1/C  FOR  LARGE  C  (08) 

IF  <  C . GT •  8 «  DO )  GO  TO  100 

IF  ( C  .  GT  .  6  .  DO  .  AND  *  M  ♦  LT .  4 )  GO  TO  200 

COMPUTE  COEFFICIENTS  FOR  CSG  EXPANSION 

IF ( L . GT . ( M+ 1 ) .  AND . C . GT .5. DO)  GO  TO  30 
LAMi-RL*  ( RL  +  i  »  E*0  ) 

LAM2=.5*(1 .-(RM2-1 . )*(RM2+i . )/( (RL2-1. >*(RL2+3. ) ) ) 
LAM3=.5*(RL-RM-1 . ) * ( RL-RM ) * ( RL+RM- 1 . >*<RL+RM>/ 

1( (RL2-3. ) * ( RL2-1 . )*<RL2-1 . )*(RL2-1. )*<RL2+1 ,))- 
2.5#(RL-RM+1» )  *  ( RL-RM+2 . ) * ( RL+RM+1 ♦ ) # (RL+RM +2 . )/ 

3 ( ( RL2  + 1 . )*(RL2+3. )**3*(RL2+5.  ) ) 

LAM4=(4,*RM*RM-1 .  )*(  (  RL-RM  +  1  .  >*<  RL-RM  +  2  .  >*(  RL  +  RM  + 1 ,  )*(RL+RM+2.  ) 
l/( (RL2-1 . )*(RL2+1 , ) * ( RL2+3 . ) **5* ( RL2+5 . )*(RL2+7. ) )- 
2  (  RL-RM- 1 .  )*(RL-RM)*(RL+RM-1.  )  *  ( RLFRM  ) 

3/ (  (RL2-5. ) * ( RL2-3 . ) *( RL2-1 . )**5*(RL2+1 . )#(RL2+3.  )  )  ) 

CL  =  LAMH-CSG*(LAM2f-CSG*(LAM3+LAM4*CSG)  ) 

EIG2=EIG3 

£IG3=EIG4 

EIG4=EIG5 

RETURN 

0  CONTINUE 


X 
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C  EXPANSIONS  BASED  ON  PREVIOUS  EIGENVALUES 
C  FIRST  THROUGH  THIRD  ORDER 
C 

IF( L . GT . ( 3+M ) )  GO  TO  50 
IF ( L . GT . ( 2+M > )  GO  TO  40 

FIRST  ORDER 

CL=2.D0*EIG5-£IG4 
1  EIG3=EIG4 

EIG4=EIG5 
RETURN 

SECOND  ORDER 

0  CL=3.D0*(EIG5-EI64)+EIG3 

1  EIG2=EIG3 

GO  TO  31 

THIRD  ORDER 

0  CL=4.D0*<EIG3+SIG5>-6.B0*EIG4-EIG2 

GO  TO  41 
00  CONTINUE 

IC=C 

IF  M>6  THEN  THE  EIGENVALUES  ARE  VERY  REGULARLY  SPACED 

IF<M.GT.6.AND.L.GT.M+1)  GO  TO  30 

USE  L*(L+1)  APPROXIMATIONS  FOR  L=M  AND  L=M+1 

IF(M.LT.  <10  +  10  )  GO  TO  200 
CL=RL*  CRL+1 . DO ) 

EIG4=EIG5 

RETURN 

COMPUTE  ESTIMATE  WITH  ASYMPTOTIC  EXPANSION 

00  B0SH=1 » DO 

IF(M.EQ.O)  B0SH=0 . DO 

CL=<R+R+1 . )*C-.25*<2.*R*R+R+R+3 . )-<R+R+l . >*<R*R+3. )/(C*16. > 
1+RM*RM+(RM-1 , >*<RL-RM)*BOSH 
EIG2=EIG3 
EIG3=EIG4 
EIG4=EIG5 
RETURN 
END 
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SUBROUTINE  CONVER < L > M » C ? CL > EIG3 > EIG5 ) 

DOUBLE  PRECISION  BLIST 

DOUBLE  PRECISION  CL»CLL»  CLSPAC>  CLU  t CORA  t CORB 

DOUBLE  PRECISION  DE,DEC>DL 

DOUBLE  PRECISION  EIG3 > EIG5 > ENR t ENRC 

DOUBLE  PRECISION  FL 

DOUBLE  PRECISION  GL > GLIST 

DOUBLE  PRECISION  C»CSQ 

DOUBLE  PRECISION  RM>RM2>R 

DIMENSION  BLIST<  250 ) » GLIST ( 250 ) ,ENR(250> 

COMMON  /BLK2/  BLIST » GLI ST » ENR 

COMMON  /BLK3/  DEC»  NDEC 

CALCULATE  SOME  HEAVILY  USED  CONSTANTS 


CSQ=C*C 

RM=M 

RM2=2.D0*RM 

DETERMINE  EIGENVALUE  SPACING  OF  PREVIOUS  EIGENVALUES 


GL-EIG5 

IF(L.EQ.M)  CLSPAC=CL 

IF < L . NE . M )  CLSPAC^E IG5-EIG3 

IU6=(L-M)/2 

IF  EST<LAST  EIG  THEN  EST =LAST  EIG 

IF(CL.LT,GL)CL=GL 

FL=CL 

JNDE=0 

IX=L-M-2*IU6 

ISC=2+IX 

LIM1=2*<L-M+C+NDEC) 

J=1 

IF(LIMl.LT.ISC)  GO  TO  25 

COMPUTE  BETA  COEF, 

DO  20  I  =  ISC  >  LIMl r  2 
R=  I 

BL I  ST ( J ) =  R*  < R- 1 , DO ) *  < RM2+R ) * ( RM2+R-1 . DO ) 

1  *CS0*CSQ/< ( RM2+2 . DOSR-1 ♦ DO ) 

2  * ( RM2+2 , D04R-1 . DO ) * < RM2+2 . D0*R-3 . DO ) * ( RM2+2 , D0*R+1 . DO ) ) 
J  =  J  +  1 

0  CONTINUE 

5  J=1 

ID21=ISC-1 

LIM11=LIM1H 

IF ( LIM1 1 . LT . ID21 )  GO  TO  35 
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C  COMPUTE  THE  GAMMA  COEF. 

C 

DO  30  I=ID21 » LIM11 1 2 
R  =  I  - 1 

GLIST( J)=<RM+R)*<RM+R+1 . DO ) + . 5D0*CSG* < 1 .DO 

1  -  ( 4  •  DOtRM-tRM-1  ♦  DO )  /  ( ( RM2  + 

2  2 . DOfcR-l •  DO ) $ ( 2 . D0&RM+2 . D0#R+3 . DO ) ) > 

J  =  J  +  1 

30  CONTINUE 

35  IFC=1 

IBLIM=LIMl/2-IX 

IGLIM=IBLIM+1 

IRI0=IU6+1 

IWl=IU6+2 

40  ENR(1)=CL-GLIST<1> 

IF  <  IIJ6  .  LT  *  1 )  GO  TO  55 

C 

C  EVALUATE  THE  CONTINUED  FRACTION 
C 

DO  50  1=1 » IU6 

ENR(I+1)=-BLIST<I)/ENR(I)-GLIST(I+1)+CL 
50  CONTINUE 

55  ENR(IBLIM)=-BLIST(IBLIM)/(GLIST(IGLIM)-CL) 

IU15*IBLIM-1 
IP=IW1+IW15 

IF  < IU15 . LT . IU1 )  GO  TO  65 
C 

C  EVALUATE  THE  CONTINUED  FRACTION 
C 

DO  60  I=IUlfIW15 
IPI=IP-I 

ENR(IPI >=-BLIST(IPI)/<GLI8T< IPI  +  1 )-CL+ENR< IPI  +  1 ) ) 
60  CONTINUE 

65  ENRC=-BLIST(IRI0)/(GLIST(IRI0+1)-CL+ENR( IRIO+1) ) 

DE=ENRC*ENRC/BLIST< IRIO) 

CORB=DE 

IF ( IBLIM . LT . I W1 )  GO  TO  75 
C 

C  COMPUTE  THE  DENOMINATOR  IN  THE  BOUUKAMP 
C 

DU  70  I=IU1 > IBLIM 

DE=ENR ( I ) SENR ( I ) /BL 1ST ( I ) *DE 
CORB=CORB+DE 

IF  ( DABS ( DE/C0RB )  .  LT . DEC )  GO  TO  75 
70  CONTINUE 

75  CORA= 1 . DO 

BE=i . DO 

IFUU6.LT. 1)  GO  TO  90 
C 
* 
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C  COMPUTE  THE  DENOMINATOR  IN  THE  BOUWKAMP 
C 

DO  80  1  =  1, IU6 

DE=BLIST( IRIO-I ) / ( ENR ( IRIQ-I ) *ENR ( IRIO-I ) ) *0E 
CQRA=CORh+DE 

IF ( DABS ( DE/CORA ) » LT . DEC  )  GO  TO  90 
80  CONTINUE 

C 

C  COMPUTE  THE  CORRECTION  TO  THE  EIGENVALUE 
C 

90  DL= ( ENRC-ENR ( IRIO  > ) / ( CORA+CORB ) 

CL=CL+DL 

C 

C  EIGENVALUE  ACCURATE  ENOUGH? 

C 

IF ( DABS  <  DL/CL ) .LT .DEC)  GO  TO  100 
IFC=IFC+ 1 

IF ( IFC*LT  »NDEC)  GO  TO  40 
100  CONTINUE 
C 

C  IF  RANGE  OF  EIGENVALUE  ALREADY  ESTABLISHED  ( JNDE  NE  0)  BRANCH 
C  TO  TEST  WHETHER  THE  RESULTING  EIGENVALUE  IS  WITHIN  THE 
C  ESTABLISHED  RANGE 
C 

IF ( JNDE  ♦  NE  *  0 )  GO  TO  120 
C 

C  TEST  OF  THE  EIGENVALUE  FOR  JNDE  =  0 
C 

IF(CL.GT.GL)  GO  TO  115 
C 

C  CONVERGED  TO  THE  NEXT  LOWER  EIGENVALUE  OF  THE  SAME  PARITY 
C 

CLL=FL 

CLU=FL+CLSPAC*1 ,5D0 
GO  TO  130 

115  IF< (CL-FL) .LT. <FL-GL) >  GO  TO  140 
C 

C  CONVERGED  TO  THE  NEXT  HIGHER  EIGENVALUE  OF  THE  SAME  PARITY 
C 

CLU=FL 

CLL=.5D0*(FL+GL) 

GO  TO  130 
C 

C  NARROWING  OF  RANGE  IF  THE  CORRECT  EIGENVALUE  HAS  NOT  BEEN  OBTAINED 
C 


120 

IF(CL.GT.CLL) 

GO 

TO 

1  no 

CLL=FL 

GO  TO  130 

122 

IF(CL.LT.CLU) 

GO 

TO 

140 

C  L  U  =  F  L 

* 
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EIGENVALUE  IS  NOW  SOMEWHERE  IN  THE  RANGE  ESTABLISHED  ABOVE 
CHOOSE  THE  MIDPOINT-  AND  REPEAT  THE  BOUUKAMP  PROCEDURE 

30  CL=  *  5D0#<  CLL+CLU ) 

FL=CL 
IFC=  1 

JNDE”JNDE+'l 

IF  MORE  THAN  50  NARROWINGS  ARE  REQUIRED  THEN  SOMETHING  IS  WRONG 

IF ( JNDE ♦ EQ  *  50 )  GO  TO  900 
GO  TO  40 
40  EIG5=CL 
RETURN 


ERROR  PRINTOUT 


900  PRINT  999, L 

999  FORMAT ( IX, 37HERR0R  IN  EIGENVALUE  ROUTINE  CONVER  AT/ 

113H  EIGENVALUE  *,I5,29H  THIS  VALUE  MAY  BE  INACCURATE) 
RETURN 
END 


38 


SUBROUTINE  OUTPUT ( A > IB , B » NACC , NUM ) 

DIMENSION  A<3) »IB(3) »B(3) »DIB(3) »DB(3) * IIB(3) t IP 1(3) »IP2(3) t 
1IP3(3) » ISIG (3) > NACC (3) 

INTEGER  PLUS  j  MINUS 
DATA  PLUS/ 1 H+/ » MINUS/ 1H-/ 

DOUBLE  PRECISION  ArB»DIBrDD 
DO  50  1=1 t NUM 

2  IF(B( I > .LE.-999.D0)  GO  TO  100 

3  IF ( B ( I ) .GE.999.D0)  GO  TO  200 
DIB(I)=IB(I> 

I  IB ( I )  =B ( I ) 

DB(I)=IIB(I) 

B  < I )  =B  < I ) -DB ( I ) 

B<I>=10.D0**B<I) 

IF < B < I ) .  GT  ♦ .99999999 DO)  GO  TO  10 
BU)=B<I)*10.D0 
I IB  < I )  =  1 IB ( I ) -1 
GO  TO  20 

10  IF ( B ( I ) .LT. 9. 9999999 DO)  GO  TO  20 

B(I )=B<I)/10.D0 


IIB(I)=IIB(I)+1 

20  B(I)=B(I)*DIB<I) 

ISIG  < I ) =PLUS 

IF ( IIB ( I ) *  GE  *  0  )  GO  TO  30 
I  IB  <  I )  »-l>XI  IB  <  I  > 

ISIG ( I ) =MINUS 

30  IP1<I)=IIB<I)/100 

IP2<I)=IIB(I>/10-IP1<I)*10 

IP3(I)=IIB<I)-IP1(I>*100-IP2(I)*10 

50  CONTINUE 

PRINT  1><A<I)»B<I)»ISIG(I)»IP1(I)»IP2<I)»IP3<I) t NACC ( I ) » 1  =  1 >  NUM ) 
RETURN 

100  B  < I ) =  0  *  D 0 

ISIG(I )=PLUS 
IP1(I)=0 
IP2 ( I ) =0 
IP3 ( I ) =0 
GO  TO  30 

200  B ( I ) =9  » 99999V9D0 

ISIG  (  I )  =PLL'S 
IP1 < I ) =9 
IP2 ( I ) =9 
IF'3  ( I  )  =9 
NACC ( I ) =0 
GO  TO  50 

1  F0RhAT(3(5XrF8'3f lXf F 10.7* lHUfAlf 311 f3Xf lit 4X)  ) 

END 

* 


39 


